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Abstract 

We present the results of a combined metadynamics-umbrella sampling investigation of the 
puckered conformers of pyranoses described using the GROMOS 45a4 force field. The free energy 
landscape of Cremer-Pople puckering coordinates has been calculated for the whole series of a 
and (3 aldohexoses. We show that the 45a4 force field parameters fail in reproducing proper free 
energy differences between chair conformers for many of the inspected monosaccharides. In the 
extreme cases of galactose, mannose and allose, the experimentally non-detectable inverted chair 
conformers become even substantially populated. The opposite behavior is observed in the case of 
idose, which is the only experimentally known aldohexose that shows equilibrium between chair 
and inverted chair. We suggest a modification to the GROMOS 45a4 parameter set, which improves 
considerably the accordance of simulation results with theoretical and experimental estimates of 
puckering free energies. 
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I. INTRODUCTION 



Within the framework of classical force fields, the number of computer experiments on sac- 
charides has grown considerably in recent years, and various systems have been addressec^i"— . 
Devising a realistic model of monosaccharides is obviously a decisive step in order for carbo- 
hydrates simulations to have enough predictive power. The accurate description of monosac- 
charides with classical force fields is not an easy task, because of the delicate interplay of 
different factors such as the presence of a high number of intramolecular hydrogen bonds, 
the competition of these hydrogen bonds with water-sugar ones and important steric and 
electrostatic effects between ring substituents in spatial proximity (see for example Ref. |28 



and references within). The problem of reproducing some carbohydrates peculiarities, such 
as the rotameric distribution of the hydroxymethyl group or the anomeric and exo-anomeric 
effects, have been addressed in various force fields, and the reader can find some comparative 



analyses in Refs. 



29-31 



However, considerably less attention has been devoted so far to the 
correct reproduction of the ring conformational properties. 

Cyclic monosaccharides can appear indeed as puckered rings, and their conformational 
transitions dramatically alter the equilibrium properties of both single sugar rings as well as 
those in oligo and polysaccharides^^. Despite the large number of possible puckered conform- 
ers, many biologically relevant monosaccharides in the pyranoid form appear almost always 
in one stable puckered conformation, the second most populated state being so unlikely not 
to be detectable by actual experimental techniques. In spite of that, several authors reported 
an inappropriate percentage of secondary puckered conformations^"— when modeling car- 
bohydrates using classical force fields such as GROMOS^^'^'^*"— or opls-aa^^ for simulations 
in solution. Regarding the latest GROMOS parameter set for carbohydrates (45a4)^^ (which 
incidentally was also tuned to improve the description of puckered conformers with respect 
to older versions), non-chair conformers have been shown to be accessible during equilibrium 
simulation runs of /3-D-glucose, with a statistical frequency of 0.7%^^. A more quantitative 
analysis has been performed recently, by computing the free energy difference between ring 
conformers of /3-D-glucose^^. The prediction of the inverted chair free energy, 16.5 kJ/mol, 
is almost 10 kJ/mol lower than most theoretical and ab-initio simulation results. 

A reparametrization of those force fields which present unphysical ring conformers at 
equilibrium would be therefore highly desirable. A quantitatively accurate description of 
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puckering free energies is a key requirement to perform proper simulations of those sugars 
which show an equihbrium between the two chair forms, but also for out-of-equilibrium sim- 
ulations of any sugar. Simulated pulling of polysaccharides, for example, has been proven to 
be an important tool for the interpretation of single molecule force-spectroscopy datai^iiMil^^ 
although the ring conformational transitions simulated using three different force fields 
(AMBER94^,AMBER-GLYCAMii and CHARMM-Parm22/SU01-^ ) led to different in- 
terpretation of the same experimental data^i. Obviously, having control on the reliability 
of force fields in reproducing puckering properties would certainly improve all these kind of 
approaches. 

Even though the importance of reproducing correctly the population of ring conformers 
is quite broadly recognized, few attempts were made so far to investigate the problem in a 
systematic way and to understand the origin of the discrepancy between simulation results 
and experimental and theoretical findings (see for example Refs. 3^,35 49). For this reason, 
we decided to study, employing a combined metadynamics - umbrella sampling approach^, 
the puckering properties of simple a and /3-D-aldopyranoses. We modeled the saccharides 
using he 45a4 carbohydrate parameter set^ of the GROMOS force field, which shows incorrect 
population of conformers for many sugars, and we identified a new set of parameters that 
satisfactorily reproduces the free energies of the main ring conformers. 

This paper is organized as follows: in Sec. [TTl the problems related to the determi- 
nation of puckering free energy landscape will be presented; in Sec. IIIII the combined 
metadynamics-umbrella sampling will be discussed, as well as the use of Cremer-Pople 
puckering coordinates^*^, their peculiarities as collective variables, and the simulation de- 
tails; in Sec. llVl the puckering free energy profile obtained with the GROMOS 45a4 force field 
for /3-D-glucose is discussed in detail and compared to recent similar calculations^'^. In SeclVl 
we present the puckering free energy profiles of the whole a and /3 series of D-aldohexoses 
modeled with the GROMOS 45a4 force field, showing that the puckering free energies are so 
poorly reproduced, that the inverted chair conformers of galactose, mannose and a-D-Glc 
become significantly populated. Subsequently, we propose a change to the 45a4 parameter 
set, showing how the free energies obtained with the use of the new parameter set compare 
favorably with available theoretical and experimental data. Eventually, some concluding 
remarks are presented in Sec. IVI[ 
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II. THE PROBLEM OF PYRANOSES PUCKERING FREE ENERGY 



In this work we are focusing on a particular class of carbohydrates, namely, aldohexoses, 
saccharides of composition C6H12O6. Aldohexoses can appear in nature in the form of six- 
membered, cyclic rings, conventionally known as the pyranose form. Pyranoses are charac- 
terized by the presence of five chiral centers located at the five ring carbon atoms. Following 
the convention of Fig. [H these chiral carbons are here denoted as CI (the anomeric carbon 
atom), C2, C3, C4 and C5 (the configurational carbon atom). This leads to quite a high 
number (2^ = 32) of diastereoisomers, characterized by the axial or equatorial orientation 
of ring substitutents. Based on the chirality at CI and C5 it is possible to classify pyranoses 
into a- and /3-pyranoses or L- and D-pyranoses, respectively^. For example, a-D-Glc and 
/3-L-Glc differ from /3-D-Glc only by the chirality at CI and C5, respectively. 

For each of these two chiralities, 2^ = 8 stereoisomers remain, which are conventionally 
named^ glucose (Glc), galactose (Gal), mannose (Man), allose (All), altrose (Alt), talose 
(Tal), gulose (Gul) and idose (Ido). The position of substituents at the chiral centres are 
reported for convenience in Tab. HI for /3-D-aldopyranoses. 

Each of these stereoisomers, in turn, can adopt different ring conformations. According 
to lUPAC recommendations^, the conformation with two parallel ring sides (four coplanar 
atoms) and the other two atoms at opposite side of the ring plane is called chair. Chair 
conformers occur in two forms: the one with C4 and CI respectively above and below the 
plane defined by the other four atoms, is identified by the symbol ^Ci, while the conformer 
showing CI above the ring plane and C4 below it is denoted by the symbol ^C4 and often 
called inverted chair. Other puckered conformers with four coplanar atoms exist, such as 
the so-called hoot and skew-boat (the latter also known as twist-boat). In the case of boats, 
two parallel ring sides define the ring plane and the atoms at the extrema of the ring are 
either above or below the ring plane (indicated for examples as '^'^B or B3 0); in the case of 
skew-boats the ring plane is defined by three adjacent atoms and the nonadjacent one with 
the other ring atoms at opposite side (e.g. ^Si) of the ring plane. Some ring conformations 
are depicted in Fig. |2J 

Conformers other than the one listed so far, which are of some relevance, are also given a 
conventional name, such as in the case of envelopes or of half chairs, and the reader can find 
a exhaustive list in Ref. 154- Obviously, many of the possible conformations are not likely to 
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contribute significantly to the chemistry and physics of monosaccharides. For pyranoses, only 
the ^Ci and conformations happen to be eligible to be the most stable conformations, 
that is, global minima which are well separated in energy from neighbouring local ones. 
Other conformers such as boats and skew-boats can in principle appear as local, scarcely 
populated minima, while half-chairs and envelopes usually correspond to maxima or saddle 
points and are typical transition states between boats and chairs. 

The concept of puckering in rings (specifically, in six-membered rings) dates back to the 
studies of Sachse^i^ on cyclohexane in the late nineteenth century, and has been later ap- 
plied also to saccharidic rings^. Although it was clear that the most stable conformers could 
have been only either of the two chairs, lack of quantitative information on their relative 
abundance stimulated many theoretical approaches (for an historical review the reader can 
look at Ref. i58i ). Semiempirical approaches accounting for the Hassel-Ottar effect (insta- 
bility of 1,3-diaxial hydroxymethyl - hydroxyl conformation) and other instability factors 
introduced by Reeves and Kelly^"— allowed to understand that for most pyranoses the 
conformation, rather that the "^Ci, is the most stable one, while in few cases an equilibrium 
between the two is possible. A quantitative estimate of the free energy difference between 
puckered conformers in aqueous solution (and not, as it is sometimes erroneously reported, 
in vacuum) came first with the analysis performed by Angyal who assigned, employing again 
an empirical scheme, specific free energy contributions to the different destabilizing interac- 
tions. The various terms were derived from measured equilibria in solution of cyclitols with 
their borate complexes^. Later, employing a molecular- mechanical approach, Vijayalakshmi 
and Rao^ obtained other estimates, which were anyway compatible with the inverted chair 
free energies predicted by Angyal. 

Semiempirical methods have been in general a quite successful approach for determining 
the stereochemical properties of numerous cyclic compounds Unfortunately, experimental 
estimates of the ring conformational free energy are not feasible for many pyranoses, mostly 
because the populations of the less stable conformers are are usually too tiny to be detected 
with probes such as NMR— For some pyranoses, however, the destabilizing effects in the 
^Ci and conformations are similar and, as in the case of idose, the relative populat ion of 
these two conformers becomes an experimentally accessible quantity. Indeed, in Refs. 



62 



63 



the free energy differences estimated for a-D-Ido and /3-D-Ido were proven to be compatible 
with the experimental findings obtained by Snyder and Serianni years later— (see Tab. HB . 
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The accordance with experiment is within ~ 4kJ/mol, but these results are of particular 
significance because they correctly predicted the preference of a-D-Ido for the inverted chair 
conformer. 

The only other pyranose for which experimental measurement of inverted chair popu- 
lation via NMR could be feasible is altrose. Some NMR conformational studies regarding 
a-D-Alt in macro cycles^*^ in solution show indeed that ^Ci/C4, and ^82 have roughly the 
same population. To the best of our knowledge, however, no results about altrose monosac- 
charide conformer populations have ever been published (NMR measurements on altrose are 
currently being carried on in our laboratory, and preliminary results^ seems to be compat- 



ible with the theoretical predictions of Refs. |62| and 1631) . For the other pyranoses considered 



in this work, only lower bounds can in principle be determined from the sensitivity of exper- 
imental probes. In this sense, semi-empirical methods predict correctly the extremely high 
free energy of inverted chairs of glucose, galactose and mannose. 

Some more recent experiments involving atomic force microscope (AFM) spectroscopy 
have allowed researchers to estimate the puckering free energy of conformers different from 
the chair ones. Marszalek, Lee, and coworkers^i^i^, for example, estimated the puckering 
free energy of glucose twisted boats employing AFM pulling on dextran, an )-6) linked 
glucan. The conformational changes involved in the elongation of this polysaccharide are 
quite complex, involving ring deformation due to its intrinsic elasticity, rotations around 
the C5-C6 bond and chair to twisted boat transitions. In order to isolate these different 
contributions, the authors performed AFM pulling also on /3-(1^4)-D-Glc (cellulose) and on 
a /3-(l— >6) linked glucan, namely, pustulan. These polysaccharides are both homopolymers 
of glucose like dextran but in the elongation process, due to the particular linkages, cellulose 
presents (apart from chain flexibility) only ring deformation, while pustulan presents both 
ring deformation and rotation around the C5-C6 bond, but no transition to the twisted boat 
conformer. By subtracting the free energy differences related to the various processes, the 
puckering free energy of the glucose twist boat was estimated to be about 25 kJ/moF^. It 
is worth mentioning, that AFM spectroscopy was flrst employed to estimate puckering free 
energies of glucose boats conformers on carboxymethyl amylose (CMA) by Marszalek and 
coworkers—, and Li and coworkers^, independently. The two groups reported values of the 
boat free energy which were in agreement and in the range 15 — 18 kJ/mol. These results 
differ considerably from the estimate of twisted boat free energy obtained from the analysis 
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of force-extension curves of dextran. However, CMA — oppositely to dextran, which has no 
intrinsic torsion — adopts preferentially a pseudo-helical conformation, and this has already 
been pointed out in Ref. O to be a possible source of error when using freely jointed chain 
models to interpret data. More recently, Kuttel and Naidoc^^ suggested that in amylose 
stretching the elongation might be due to a complex rotation of glycosidic linkages, and 
that chair to boat transitions could play little role in the process. Finally we would like to 
mention that ^Ci to ^C^ transitions could be in principle investigated by means of AFM 
spectroscopy of, e.g., pectin^l. Their more complex two-step transition prevented so far an 
estimate of the puckering free energy, but it still represents — to the best of our knowledge 
— the only approach which could provide direct estimate of inverted chair free energies of 
the pyranoid form of aldohexoses other than idose and altrose. 

In the field of computer simulations, ab-initio predictions of the puckering free energy 
of pyranoses are even more scarce: we are only aware of a Car-Parrinello metadynamics of 
glucose in vacuunt^ (where, unfortunately, use of non-optimal collective variables has been 
made of^), a Hartree-Fock/coupled clusters calculation (MP2/cc-pVTZ) with inclusion of 
vibrational, rotational and solvation free energies contributions^, and a calculation at the 
B3LYP/6 — 311 + +G** leveF^. The calculations of Appell and coworkers^^ predicted a free 
energy difference for the conformer of a-D-Glc and /3-D-Glc of 20.40 and 29.18 kJ/mol, 
respectively. It has to be noticed, that calculations performed by including solvation effects, 
although only at the MP2/cc-pVTZ level^ led to a free energy of /3-D-Glc conformer 
of about 57 kJ/mol. This very estimate appears to be substantially higher than all other 
ones reported so far. In absence of similar ab-initio calculations for idose, that could be 
employed as benchmark system, it is in our opinion difficult to assess the confidence level of 
this estimate. However, the possibility that the values reported in Ref. |62| and Ref. |63| could 
underestimate the free energy of inverted chairs for stereoisomers other than idose can not 
be completely ruled out. 

Given the results briefiy reviewed so far, it is obvious that clear-cut experimental reference 
values for sugar ring puckering free energies are not yet available. The theoretical estimates 



presented in Ref.|62!and in Ref.|63|are mutually consistent, in agreement with recent ab-initio 
calculations^^, and were able to predict with reasonable accuracy the population of chairs in 
a, /3-D-Ido, and possibly altrose^. The results of AFM pulling on dextran are also compatible 
with these estimates, in the sense that the free energy of boat conformers is expected to be 
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not lower than that of inverted chairs. On the other hand, some results appear to be not 
compatible with this set of estimates. Ab ab-initio calculatiort^ assigned a much higher free 
energy to glucose inverted chairs, while other data form AFM spectroscopy of CMA—^^ go 
into the opposite direction. By and large, we feel confident that the two series^*^ represent 
reasonable approximations of pyranoses free energy difference between chair conformers, and 
in this work we will use these values as a reference when testing new force field parameter 
sets. 



III. METHODS 

A. Refining Metadynamics with Umbrella Sampling 

The term metadynamics identifies a number of techniques that have been devised dur- 
ing the last decade to accelerate dynamics for systems displaying meta-stabilities. Various, 
slightly different algorithms have been proposed in literature, that exploit the same under- 
lying idea of the first metadynamics implementation^" (sometimes referred to as discrete 
metadynamics): the Lagrangian version^-, which relied on the introduction of fictitious de- 
grees of freedom associated to the reaction coordinates; the direct metadynamics^, which 
is probably the most employed variant; the well-tempered metadynamics^, which assures 
a convergent estimate of the free energy profile. All these variants share the usage of a 
time-dependent biasing potential, Uuas[s{x),t], to ease the exploration of the phase space 
along suitably chosen collective variables, s{x), in turn themselves functions of the atomic 
coordinates x. The collective variables must represent all slow degrees of freedom charac- 
terizing the system, in order for metadynamics - as well as any other free energy profile 
reconstruction method - to be meaningful^. 

The direct metadynamics approach shares many traits with other enhanced sampling 
methods, such as the Local Elevation method-- or the adaptive umbrella samplingSS (detailed 
lists of enhanced sampling techniques based on molecular dynamics can be found in Refs. |43 
and l84j) . but the fact that the free energy landscape in metadynamics is usually estimated 
as the negative of the biasing potential in a single, out-of-equilibrium sweep, makes the 
technique susceptible, at least in principle, of errors introduced by the deposition protocol. 
Distinct approaches have been devised to estimate the statistical and systematic errors in 
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metadynamics, and also to recover a truly equilibrium free energy profile^>^i2I. Among 
these, the approach proposed by Babin and coworkers^ is in our opinion one of the most 
intuitive and versatile ones, since it allows in a simple way to simultaneously estimate the 
statistical error and to eliminate systematic errors introduced by the deposition protocol. 
The basic idea is the following: a metadynamics run is performed up to the build-up time 
tb, so that the whole collective variables space has been explored, and the total potential 
energy - sum of the physical and the bias potential - at the end of the run reads 

Vix) = U{x) + Uu^[six),h] . (1) 

At this point, a molecular dynamics simulation in the potential V is performed, much in 
the spirit of umbrella sampling, whereas the biasing potential has been determined in an 
adaptive way by the metadynamics run. The dynamics is characterized by an almost diffusive 
behavior in the collective variables space, since the meta-stabilities have been removed by 
the bias potential, provided that all the states of the new potential V{x) are separated by 
energy barriers comparable or lower than the thermal energy scale k^jT. The deviation from 
a truly diffusive behavior is due by residual features of the free energy landscape which are 
originated from the statistical and systematic errors in the metadynamics run. During this 
equilibrium run, the phase space is sampled with probability density 

g-/3{F(s)+!7bi,,[5(a;),tb]} 

that can be estimated by computing the histogram pbias(s) = (t — tb)~^ Jl^^5[s — s{t')]dt' 
during the run. Eventually, the free energy profile is given, up to an additive constant, by 

F{s) = -f/bias[s(x),tfe] - /CsTlupbiasls). (3) 

In the exact expression Eq. [3l the term kj^T In Phias{s) can be regarded as a correction to what 
is the standard metadynamics estimate of the free energy landscape, F^^tais) = —Uuas{s). 
This correction term compensates for the systematic errors introduced by the deposition 
protocol, which are not completely under control in the metadynamics run. In other words, 
from a simple metadynamics run there is no way to guarantee that the term pbias(s) has 
become a constant within statistical fluctuations. The corresponding statistical error can be 
estimated from the fluctuations of puasis) using standard error analysis. Differently from 
the approach involving only a metadynamics run, one should not worry about the speed of 
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the deposition process (which involves the height, width and deposition rate of the Gaussian 
functions) as long as the subsequent equilibrium run is ergodic. This means that the Gaussian 
functions placed during the build-up phase have only to (a) be smaller than kj^T, (b) being 
of width comparable or smaller than the finest detail of the free energy landscape which is 
deeper than kj^T and (c) cover the whole conformational space of the collective variables. 

This idea of employing the biasing potential with an umbrella-like sampling has been 
applied, besides to met adynamics^, also to the local elevation method^. 

B. Puckering Coordinates 

The generalized coordinates introduced by Cremer and Pople^i^ can be used to identify 
puckered conformations of rings with an arbitrary number of members. Their definition 
makes use of the projections zj of the position vector of each ring atom onto the axis of 
the mean ring plane. In the case of six-membered rings they can be defined in terms of the 
distances Zj clS db functions of 3 parameters. 




^3 = v1E(-i)'"^^' • (6) 

This coordinate set {q2,(p2,Q3) can be conveniently expressed as a spherical coordinate set 

q2 cos (p2 = Q sin 6 cos 
\ q2 sin (f)2 = Q sin 9 sin (p (7) 
qs = Q cos 6, 

V 

where 9 E [0,7r], G [0,2n) and Q, the so-called total puckering amphtude, is defined as 

There are two main interconversion paths in pyranoses, namely, the inversion path — 
connecting the ^Ci and conformers — and the pseudo-rotation patb^ — connecting 
the more flexible boat and skew-boat conformers — that can be easily represented in this 
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coordinate set. The inversion path develops along the 9 coordinate from the ^Ci conformer 
at = to the ^C^ one at = tt, while the pseudo-rotation one develops along 0, at ^ = 7r/2. 
Notice that is a 2-7?— periodic coordinate, meaning that points at = or = 27r (at a 
given value of 9) represent the same conformer. This is not true for 9, which is not periodic. 

In this sense, the spherical set (Q, 9, 0) has the advantage that only the two angular 
variables are needed as collective variables in order to perform a proper - that is, ergodic 
and unbiased - exploration of the puckered conformations space of typical six-membered 
rings. This is because along the radial direction no meta-stabilities occur, and the radial 
coordinate Q relaxes fast enough to be ergodic for every reasonable set of potential (that 
is, when the bond lengths are rigid or quasi-rigid). As we showed in a previous work^^, not 
every representation of the Cremer-Pople coordinates is equivalent to the end of being used 
as collective variables for a conformational search. In particular, any two-dimensional subset 
of Cremer-Pople coordinates whose functional form involves also biasing forces along the 
direction of Q might suffer from lack of ergodicity and, therefore, lead to biased sampling^. 
In this work we make only use of the angular variables of the spherical coordinate set as 
collective variables for the metadynamics run. 

Alternative generalized coordinates can be in principle employed to characterize the puck- 
ered conformers of six-membered rings, such as the three out-of-plane dihedrals introduced 
by Strauss and Pickett^, or other definitions based on three internal dihedral angles^"—. All 
these alternative schemes produce good puckering coordinates, in the sense that they allow 
to uniquely map the complete puckering conformational space. However, in the context of 
conformational search using accelerated methods, oppositely to Cremer-Pople coordinates, 
(a) they require the exploration of the complete three-dimensional space, being therefore less 
convenient than a two dimensional phase space search, and (b) the set of micro-states cor- 
responding to the different conformers generally define less simple surfaces which, as it will 
be discussed later, might affect the determination of the relative populations of conformers. 



C. Simulation Details 



The metadynamics and equilibrium simulations have been performed using a version of 
the GROMETA simulation package^"—, previously modified to allow the usage of the 9 and 



angular coordinates of Cremer and Pople, as described in Ref. 



77. For each of the 16 
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D-aldopyranoses, a system composed of the respective sugar ring in a cubic simulation box 
filled with 504 water molecule was set up. The SPG^Si model has been used to describe 
the water molecules employed to hydrate each of the pyranoses. The settle algorithm has 
been used to make water rigid^^, and all bond lengths in the sugar molecules have been 
constrained using the shake algorithm^*^. An integration step of 0.2 fs has been used for 
every phase of the simulations. Before starting each run, a 100 ps long molecular dynamics 
simulation with no bias has been performed to equilibrate the different sugars in their ^Ci 
conformer. 

The metadynamics part of the run consisted in a 4 ns long simulation, using Gaussian 
potential functions of height 0.5 kJ/mol and with 0.05 rad for both angular variables, to 
build up the biasing potential. Gaussian functions have been placed every 200 integration 
steps. The Nose-Hoover— thermostat and the Parrinello-Rahmani^i barostat have been 
applied to simulate isothermal-isobaric conditions at 300 K and 1 atm using relaxation times 
of 0.1 and 1.0 ps, respectively. The simulation box edges have been kept orthogonal, and 
have been rescaled using an isotropic pressure coupling, which controlled the trace of the 
pressure tensor. 

Each metadynamics run has been then followed by a 4 ns long equilibrium molecular dy- 
namics run at the same thermodynamic conditions, using the set of Gaussians placed during 
the metadynamics to generate the time-independent biasing potential. At a difference with 
Ref. |50|, where the Lagrangian metadynamics with truncated Gaussians has been employed, 
we make use of the standard direct metadynamics. During this run, the histogram pbias(6', 0) 
has been collected by samphng configurations every 40 fs on a grid of 60 x 60 points. 



IV. PUCKERING FREE ENERGY OF GLUCOSE 

The combined metadynamics-umbrella sampling has been first applied to the calculation 
of the puckering free energy profile of /3-D-Glc, employing the standard GROMOS 45a4 pa- 
rameter set^. In Fig. [3] we present the free energy profile as a function of the Cremer-Pople 
angular variables 6 and 0, showing isolines separated by 10 kJ/mol (upper panel). In order 
to facilitate the comprehension of the plot, the projection of the free energy profile onto the 
= plane (lower panel) is also provided: the lower contour of the colored region in the 
projected profile allows to easily identify the minima along the 9 coordinate (corresponding 
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to the minima of this lower contour) and the transition states, or free energy saddle points 
(corresponding to the maxima of the lower contour). On the (6', 0) plane, ^Ci (0,-), ^C^ (vr,-), 
and ^'°B (7r/2,0) conformers are clearly recognizable as minima basins. Notice that, due to 
the periodic nature of (p, the '^'^B conformer appears in this representation to be split in two 
across the = line. Another local minimum basin, more shallow than the previous ones 
is located around the ^83 (7r/2,77r/6) conformer and seems to include also some other near 
boat-like conformer. This kind of occurrence can be a natural feature for gluco-pyranoses, 
because of the high flexibility of the ring conformers along the pseudo-rotational path (lo- 
cated aX 6 = 7i/2). A more detailed description of this local minimum is anyway beyond the 
interest of this work, because of its very high free energy. 

In addition to the free energy profile, in Fig. H] we show an analogous plot, presented 
both on the {9,(f)) plane (upper panel), and projected onto the = plane (lower panel), 
of the correction to the free energy profile coming from the umbrella sampling phase, 
— /cgT Inpbias(s). Performing a screening of this contribution is quite important, because 
it allows to check if residual meta-stabilities (recognizable as basins) are left after the meta- 
dynamics phase, and whether an ergodic sample of the interesting region in the {6, 0) plane 
has been performed or not: poorly sampled regions will appear wither as peaks. In Fig. HJit 
is possible to see that the metadynamics run performed reasonably well. Indeed, the paths 
connecting the original metastable states do not show residual metastabilities, being the 
more pronounced minimum, located close to (0,7r/3), only ^ 2 kJ/mol deep. It is seen that 
the complete space of puckering angles has been sampled properly. Corrections to the free 
energy difference between ^Ci and are ~ 2 kJ/mol. In case of the next most populated 
state, '^'^B, the correction increases to ~ 4 kJ/mol, while for other states the correction can 
reach, but never exceeds, ~ 6 kJ/mol with respect to the free energy of ^Ci. The umbrella 
phase thus contributes in a significant way to the estimate of free energy differences, and 
— as expected — appears to be more important for less populated states and for transition 
states. 

Both the and ^'*^B puckering free energies (local minima at about 10 and 17 kJ/mol, 
respectively) appear to be lower than what one would expect. In particular, the free energy 
of the inverted chair is about 15 kJ/mol lower than both our reference values^i^. This free 
energy difference leads to a inverted chair population of about 1%, which can be observed 
during regular molecular dynamics runs at equilibrium. While not strictly incompatible with 
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NMR experiments, which usually have a resolution of ~2%, this value is certainly strikingly 
lower than all other theoretical estimates. This fact is even more important if one considers 
that /3-D-Glc is supposed to have the largest free energy among aldopyranoses. It is 
thus expected that the conformers of the other stereoisomers could be characterized 
by even lower free energy differences. As it will be discussed in SecJVl this scenario is only 
partially correct, as unexpected patterns of the chair-inverted chair free energy difference 
appear along the series. 

While the value of the free energy at the local minima is certainly interesting, a more 
direct connection with experimentally measurable quantities is performed through the cal- 
culation of the population of the basins associated with each of the recognizable conformers. 
To this aim, the {6, 0) plane has been partitioned in four regions covering the four main 
basins: 9 G [0,7r/3], associated to ''Ci; 9 G [27r/3,7r], associated to ""^04; 9 G [7r/3,27r/3] and 
G [— 27r/3, 27r/3], associated to ^'°B; 9 G [7r/3,27r/3] and G [27r/3, 47r/3], associated to 
the mixture around ^83. This partitioning of the conformational space appears to be quite 
natural, since the separating lines are located with good approximation along the maxima 
of the free energy surface (as is evident especially from the side-view of the free energy 
landscape, lower panel of Fig. |3]) and, therefore, also close to the transition states. The 
population of a region S can then be calculated simply as 



is extended to the whole [9, 0) space. According to this numerical estimate, the ^Ci,^C4,^'*^B 
and ^83 conformations account for 98.42, 1.52, 0.0461 and 3.7x 10^"^% of the total population, 
respectively. 

When the manuscript was in preparation, we discovered that Hansen and Hiineneberger 
already performed an analysis of the puckering free energy of glucose^. To our surprise, 
their estimates of the inverted chair free energy differ by a non-negligible amount from 
that presented here, although the same force field and thermodynamic conditions have been 
used. It has to be noticed that in the approach of Hansen and Hiinenberger, free energies are 
derived from populations by inverting Boltzmann formula. We did the same and obtained a 




(8) 



where the integral of the normalization factor 




(9) 
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value of 10.43 kJ/mol for the free energy of ^C^, while in Ref. |43| this value ranged from 16.0 
to 16.5 kJ/mol (see Tab. 4 in Ref.—) The difference between these two approaches is quite 
important, as it translates to inverted chair populations of ~ 1% and ~ 0.1% for the results 
of this work and of Ref 4^, respectively. 

We give here a tentative explanation of this difference. It has to be noticed that the 
difference is most probably not related to which set of collective variables (Pickett dihedrals 
versus Cremer-Pople coordinates) has been used to enhance the sampling, because both 
represent every slow degree of freedom and allow to span the puckering conformational space 
in an ergodic way. Systematic errors related to the accelerated dynamics should have been 
eliminated in both approaches, since both sampling methods provide an unbiased estimate, 
thanks to the equilibrium sampling. 

Among the possible reasons are the use of different algorithms for the simulation of 
isothermal-isobaric conditions, the different system size, cut-off and long range corrections. 
Another explanation could be related to the way different states are defined: partitioning the 
conformational space using the three Pickett angles might lead to a miscounting of states, 
when using convex shapes like cubes or ellipsoids to perform cuts, as free energy basins in the 
three-dimensional space of Pickett angles can have have more complex, non-convex shape 
(compare for example Figs 2c and 3c in Ref. l43l ). By employing Cremer-Pople coordinates 
only two coordinates, 6 and (p, need to be taken into account. Performing the appropriate 
cuts to identify the different states is far a simpler task in two dimensions than in three, 
also taking into account that only straight cuts along the 6 and directions are needed to 
identify in a proper way the different basins. Still, more investigation would be needed to 
address properly this subject, which is ot of the scope of this work. 



V. PUCKERING FREE ENERGY OF ALDOPYRANOSES 



A. The 45a4 Parameter Set 



Apart from the problems related to a correct definition of different puckered states, it is 
clear that the calculation of the free energy profile of /3-D-Glc alone cannot be considered 
exhaustive nor satisfactorily. The theoretical estimates predict a great variety in the confor- 
mational free energies of the 16 stereoisomers of /3-D-Glc, and puckering properties have to be 



15 



checked separately for each of them. For these reasons we decided to compute in a systematic 
way the puckering free energy landscape for the whole series of a and /3-D-pyranoses. 

From the point of view of the force field, the stereoisomers of glucose differ only slightly, 
namely (a) the order of the two central atoms involved in an improper dihedral interaction 
at a chiral centre has to be inverted, in order to move a residue form the equatorial to the 
axial conformation; (b) a anomers are distinguished from /3 anomers by different torsional 
interactions on the 05-Cl-Ol-Hl dihedral, and (c) those sugars having 04 and C6 located 
on the same side of the ring plane (galactose, talose, gulose, idose) are modeled with different 
parameters for the 05-C5-C6-06 and C4-C5-C6-06 dihedral angles^. 

Simulations of the remaining 15 steroisomers of glucose have been performed using the 
same protocol employed for /3-D-Glc. We summarized the results in Fig. Oand Tab. IIIII In 
Fig. |5]we report the free energy difference between the and ^Ci conformers of a and 
/3-D-pyranoses modeled using the 45a4 force field, along with the theoretical estimates of of 
Angyal^ and of Vijayalakshmi and coworkers^. Two horizontal dashed lines are also drawn 
at and 5 kJ/mol (approximatively 2k^T at room temperature), highlighting the thresholds 
below which the inverted chair population becomes greater than the chair one, and below 
which the inverted chair population becomes noticeable, respectively. In Tab. Illlt free energy 
differences and populations for the complete a and /3 series are listed. Beside those of ^C4, 
the free energy difference and population of the next leading conformer are also reported, as 
well as the location on the {Q^ (p) plane (and the closest recognizable conformer) of the first 
transitions next to ^Ci. Concerning the population of the leading next conformer, the values 
reported were calculated on the 6 G [7r/3,27r/3] region. This choice takes into account all 
other local minima present along the equatorial line, but their free energy is always so large 
(as it has been seen for /3-D-Glc), that this approximation doesn't change substantially the 
population of the next leading conformer. 

The differences between the theoretical estimates and the simulation results obtained 
using the GROMOS 45a4 force field are striking. First of all, none of the 16 sugars investigated 
presents a chair/inverted chair free energy difference in quantitative agreement with the 
theory, as differences are usually larger than 5 kJ/mol, and in many cases even larger than 
10 kJ/mol. More importantly, many of these values are in marked qualitative disagreement 
not only with the theoretical results, but also with experimental evidence. In fact, a-D- 
Glc, a-D-Gal, a-D-Man, /3-D-Gal and /3-D-Man present an inverted chair free energy which 
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is lower than 5 kJ/mol (approximatively 2 /cg^ at room temperature). This means that a 
sensible population of inverted chairs, of the order of 10% is expected at equilibrium, in 
contrast with no experimental evidence of the occurrence of this conformer. The behavior of 
-D-Tal is even more pronounced, displaying an inverted chair free energy close to zero or, for 
the /3 anomer, even negative. This should have been the case of idose, whose inverted chair 
conformers have been experimentally detected. On the contrary, the puckering free energy of 
idose inverted chairs simulated using the 45a4 set of parameters results to be greater than 10 
kJ/mol, therefore ruling out the possibility of observing idose inverted chairs in equilibrium 
simulations. 

The GROMOS 45a4 force field appears to be unable not only to compare quantitatively 
with experimental and theoretical results, but — even more importantly — to reproduce 
the qualitative behavior of any of the two series. Given the ubiquitous presence of galactose 
and mannose in relevant oligo and polysaccharides of biological origin, the inability of the 
force field to prevent appearance of inverted chairs at room temperature seems to be a severe 
drawback, at least for out-of-equilibrium simulations. 

While the free energy of different ring conformers is certainly an important physical quan- 
tity, one should not overlook the importance of the kinetics of the conformational transitions. 
One might reason that alternate conformers might not be seen during equilibrium simula- 
tions, if the inverse transition rate is much longer than the typical time interval spanned by 
a simulation. This pragmatic approach could be hazardous, given the fast pace of increase 
in simulations sizes and lengths, but nevertheless appealing. Krautler and coworkers^^ re- 
ported that in 200 ns long simulation runs of /3-D-Glc, /3-D-Gal, /3-D-Man and /3-D-talose, 
all sugars but glucose remained for more than 99.9% of the time in the chair conformation, 
while glucose was found in boat and twisted conformation for the 0.7% of the time (giving 
a rough estimate of the characteristic time of escape from the chair conformer basin of 10 
ns). 

Although 200 ns is a time much longer than that of most simulations, one should keep in 
mind that conformational transitions are stochastic events, and a characteristic time of 10 
ns might lead to a considerable amount of "unwanted" conformers in simulation runs much 
shorter than 200 ns, but with more than just one sugar molecule in solution. We tested a 
setting which we consider to be representative of a typical simulation of medium to large 
size, namely, of a 25 ns long run at constant temperature and pressure of 512 /3-D-Glc and 
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25x10^ water molecules in a simulation box with an edge of approximately 9.6 nm, using 
a cut-off of 1.3 nm for every nonbonded interaction, and an integration timestep of 2 fs. 
Every other parameter and algorithm employed has been the same as in the metadynamics 
run discussed so far. Indeed, we observed the appearance of both boats and inverted chairs 
conformers, with a statistical frequency of 0.06% and 1.1%, respectively (see the upper panel 
of Fig. [6l These values are close to the one expected from the free energy calculations (see 
Tab. Illl|) . However, a look at the time evolution of the number of inverted chain in the 
simulation box (lower panel of Fig. [6]) tells us that equilibrium has not been reached yet. 
Still, this result is to our opinion quite valuable, because it gives some information about 
the kinetics of ring conformational transitions in /3-D-Glc, showing that it is not unlikely to 
observe inverted chair conformers in equilibrium simulations of conventional size, using the 
45a4 parameter set. 

We performed similar simulations of other sugars, namely, /3-D-Gal and a-D-Glc, and 
in both cases we observed the appearance of inverted chairs, although to a much smaller 
extent, showing that the kinetics of the chair to inverted chair transition is much slower in 
these cases (consistently with the results of Ref. l28l ) . This fact is obviously related to the 
height of the barrier that separates ^Ci conformers from other ones (see Tab. Illll) . which 
appears to be lower in /3-D-Glc than in all other cases. 



B. Force Field Re-parametrization 

After realizing the difficulties of the 45a4 parameter set in reproducing puckering prop- 
erties of pyranoses, we planned to re-parametrize the force field, by finding a minimal set of 
changes that could fix at least the qualitative aspects discussed so far. At a first glance, this 
task could seem a bit intimidating, because puckering variables (and therefore the relative 
free energy landscapes) depend directly on all six ring atoms and also — indirectly, but 
possibly to a considerable extent — on the ring substituents. The number of parameters 
on which puckering free energy depends is, therefore, quite high. A completely automated 
procedure is out of question, and one needs to adopt an heuristic approach. 

To keep the re-parametrization as general as possible, and the number of parameters to 
be tuned low, we decided that our approach should adhere to the following criteria: (a) only 
parameters not directly involved in inter-molecular interactions should be tuned; (b) the 
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changes should not be sugar- dependent; (c) the changes have to preserve previously known 
or already tuned molecular properties; (d) the inverted chair free energy of most common 
sugars (glucose, galactose, mannose) should be higher than 10 kJ/mol; (e) the trend of the 
inverted chair free energy as a function of the sugar type, as well as (f) the approximate 
offset between the inverted chair free energies of a and /3 anomers, have to be reproduced. 

Within this framework, a complete, quantitative agreement for every sugar type is most 
probably not feasible, especially given the constraints (a) and (b). However, we found that 
quite reasonable results can be indeed achieved with minimal parameter changes. Point 
(a) requires Lennard- Jones parameters and partial charges of the 45a4 parameter set to 
be preserved. Together with the requirement at point (d), this suggests that only angular 
or torsional interactions involving three or more ring atoms should be the target of our 
optimization. This is because among the known properties which are well reproduced by the 
GROMOS 45a4 force field there are the rotameric distribution of the hydroxymethyl group 
and the conformation of the glycosidic linkages (in di-saccharides). We realized quite soon 
that changing the stiffness of angular interactions did not change the puckering free energy 
noticeably. Spieser and coworkers^ showed that the height energetic barrier between ^Ci 
and can be increased by stiffening the angular interactions. However, this affects only 
the free energy of the transition states, and not the free energy difference of the metastable 
conformers. 

Therefore, we concentrated on the torsional interactions, and noticed that every torsional 
interaction involving three ring atoms (C1,C2,C3,C4,C5 or 05) and either C6 or any hy- 
droxyl oxygen (01,02,03,04) was either not present (this is the case of 04-C4-C5-05, 
C3-C4-C5-C6, 02-C2-C1-05, C1-05-C5-C6 and C5-05-C1-01) or present with a phase 
term cos(5) = +1 and multeplicity n = 2 in the torsional potential 

f/(0) = A;^ [1 + cos(5) cos(n0)] , (10) 

associated to the dihedral angle 0. As it was pointed out in Ref. |43|, such a phase favors the 
axial conformation of the substituent, with respect to the equatorial one, whereas a negative 
value of cos(5) would favor the equatorial conformation with respect to the axial one. Indeed, 
one of those interactions (02-C2-C1-05) has been eliminated in the 45a4 version of the 
GROMOS force field, for the precise purpose of stabilizing the ^Ci conformer. While this is 
true for glucose, it is not, e.g.^ for mannose, whose substituent in C2 is axial in the chair 
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conformer. Therefore, the change in the 45a4 set that stabihzed the glucose chair, acted in 
the opposite direction for mannose, stabihzing the inverted chair. 

It is clear that any change on these torsional parameters will have a profound (and sugar- 
specific) effect on the whole series of pyranoses. Therefore, to make grounded changes to the 
force field, one has first to understand how the pattern of axial and equatorial substituent 
in the ^Ci conformer, for given torsional interactions involving the five chiral centres, can 
influence the properties of the free energy curves of the a and /3 series. 

Firstly, we realized that by changing the sign of cos((5) in the C1-C2-C3-03 and C5-C4- 
C3-03 interactions, one can at the same time rise the free energy of glucose, galactose, 
mannose, and talose, and lower that of allose, gulose, altrose and idose. Looking at the 
location along the free energy series of these sugars, it appears that these interactions are 
leading candidates to recover the approximate monotonous trend in the free energy and, 
consequently, to fix consequently point (e). 

Another important role in determining the shape of the free energy curve is played by a 
dihedral interaction which is not present in the 45a4 parameter set, namely, the one involving 
the substituent at the C5 chiral center. The chirality of C5 is the same for all 16 steroisomers 
of D-Glc, and can be actually exploited to introduce a global shift for the free energies 
of the whole series of 16 D-pyranoses. While this gives some freedom in our parametrization 
process, it should be kept in mind that, most probably, parameters optimized this way are 
not valid to model the series of L-pyranoses. 

In addition to the interactions discussed so far, to reproduce correctly the gap between 
the free energy curves of the a and /3 series — as it is apparent, for example, in the 
theoretical data presented in Fig.|5] — the torsional interaction for the C3-C2-C1-01 present 
in the 45a4 set has to be modified. The chirality of the CI carbon atom differs only between 
the a and /3 anomers, and is certainly playing a role in modulating the height of the free 
energy gap between a and /3 anomers. 

The case of idose is striking, because non-bonded interactions between the substituent 
in CI and the other interaction centres are, for this sugar, relatively weak with respect to 
those of the other substituents. Therefore, the torsional interaction C3-C2-C1-01 will be 
the leading actor in determining the gap in free energy between a-D-Ido and /3-D-Ido. By 
performing a simulation with no torsional interaction on C3-C2-C1-01 for both the a and fi 
anomers of idose, the gap between the ^C4 free energies, indeed, vanished (within statistical 
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error bars of 0.4 kJ/mol). This fact could in principle be exploited to fix the C3-C2-C1-01 
interaction so that it reproduces the gap between the ^C^ free energies of a and /3 idose. In 
practice, however, we were unable to proceed along this line, because we did not manage to 
obtain satisfactorily puckering properties for idose and all other sugars at the same time. 
This could be related either to the fact that, admittedly, we did not perform a complete 
exploration of the parameter space due to the high demand of computational resources of 
this task or, to the fact that also non-bonded interactions should be re-parametrized, given 
the complex interplay typical of effective force fields. Eventually, some additional corrections 
were made by adding the same kind of torsional interaction on the dihedral C4-C3-C2-02, 
which is related to the position of the substituent at the chiral carbon C2. 

In summary, the phase and amplitude of the C3-C2-C1-01, C4-C3-C2-02, C1-C2-C3- 
03, C5-C4-C3-03 and C1-05-C5-C6 torsional interaction were tuned by trial and errors in 
order to obtain free energy curves in better agreement with experimental and theoretical 
data. The set of torsional interactions and their parameters for the proposed modification to 
the 45a4 parameter set a re reported in Tab. IIVI It is worth mentioning that the strength of 
the C1-C2-C3-03 and C5-C4-C3-03 torsional interactions had to be set to a much higher 
value (2.4 kJ/mol) than that of the other ones involving three ring atoms and one hydroxyl 
oxygen. Given the similar chemical nature of the quadruplets of atoms involved (beside that 
of the anomeric oxygen), such asymmetry appears to be peculiar. This might originate from 
other interaction terms already present in the 45a4 set, whose influence on the puckering 
free energy is not yet understood, and that might deserve a separate investigation. 

The new parametrization of the GROMOS force field for sugars was performed with the aim 
not to change properties other than puckering, and dihedral interactions directly involved in 
the rotameric distribution of the hydroxymethyl group were thus not changed. This choice 
alone, however, is no guarantee that this quantity is not affected We checked explicitly that 
these modifications did not affect the rotameric distribution of the hydroxymethyl group, 
calculating the free energy profile of the C4-C5-C6-06 torsional angle for the 45a4 and 
new sets of parameters. The results obtained with the two parameter sets did not differ 
more than 2% between each other. The free energy surfaces have been calculated employing 
the same metadynamics/umbrella sampling approach employed for the calculation of the 
puckering free energy, but using a Gaussian width of 0.1 rad. 

We summarized the results obtained using our modified set of parameters in Figs. ITllHllQ] 
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and Tab. |Vl In Fig. [7]the free energy differences between inverted chair and cliair conformers 
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and Ref. 
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Tlie improvement 



are compared again witli the theoretical predictions of Ref. 
with respect to the results obtained with the 45a4 set (Fig. E]) is striking. Both the a and /3 
series reproduce now the qualitative trend of the theoretical estimates. Galactose, mannose, 
and talose are not anymore below the 2/cbT threshold, and the value of gulose, altrose and 
idose free energies diminished considerably. Also the gap between the a and /3 anomers is now 
reproduced reasonably well, being on average ^ 5kJ/mol. Noticeably, the height of the free 
energy barrier of /3-D-Glc has been increased considerably, from ^ 25kJ/mol to ~ 45kJ/mol. 
This increase, which changes dramatically the kinetics of the conformational transition of (3- 
D-Glc, has no counterpart in any other /3— anomer. The situation for a— anomers is slightly 
different, because the height of the free energy barrier close to ^Ci increased markedly for 
a-D-Tal and a-D-Gul and, at the same time, decreased for a-D-All and a-D-Alt. 

Unfortunately, the population of inverted chairs in a-D-Ido is still lower than that of 
chairs, whereas both theory and experiment show a preference for the conformer. How- 
ever, given the fact that these changes are not sugar-specific, the result obtained is to our 
opinion still remarkable, as the ability to reproduce puckering properties has increased dra- 
matically, with respect to the 45a4 set. 

We expect that a more in-depth analysis and parameter fine-tuning could lead to even 
better agreement with theory and experiments. 



VI. CONCLUSIONS 



We performed the first systematic calculation of puckering free energy landscapes for 
the series of a and /3-aldohexoses using the latest parameter set (45a4) for the GROMOS 
force field. While the realism in reproducing the population of inverted chairs and boats for 
/3-D-glucose has surely improved with respect to the 45a3 set, it is still in marked quan- 
titative disagreement (around 10 kJ/mol) with available theoretical models. Besides that, 
we discovered that the free energy difference between inverted chairs and chairs of other 
pyranoses modeled using the 45a4 set leaves still much to be desired: galactose, mannose, 
allose, and a-D-Glc present a sizable population of the inverted chair conformer, which has 
never been observed experimentally. On the contrary, the population of the inverted chair 
conformer of idose resulted to be negligible in simulation, whereas theoretical approaches 
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and experimental results confirm the equilibrium between the two chair conformers of idose 
in solution. 

After analyzing the role that some torsional interactions have in modulating the inverted 
chair free energy curves for both the a and /3 series of D-aldopyranoses, we tested a new set of 
parameters which proved to be quite successful in recovering the trend of the inverted chair 
free energy for all 16 stereoisomers under study. Our parametrization reproduces free energy 
differences in accordance to experimental and theoretical data always within 5 kJ/mol, but 
in most cases within 2.5 kJ/mol, that is, Ik^T. A closer agreement with the theoretical model 
of Angyal^ or Vijayalakshmi and coworkers^ is probably not needed, given the uncertainties 
to which these models are subjected. Indeed, while the theoretical models do not provide 
any confidence interval, an approximate picture can be made by comparing the theoretical 
results with experimental data on idose^, which roughly fall within a 2-3 kJ/mol interval. 

The improvement attained with the introduction of this new parameter set is to our opin- 
ion substantial, and reached the goal of reproducing known puckering free energy differences 
while keeping other properties, such as inter-molecular interactions and the rotameric distri- 
bution of the hydroxymethyl group, unchanged. These modifications to the GROMOS force 
field will allow to perform more realistic simulations of D-aldopyranoses. They represent an 
certain improvement in the study of carbohydrate equilibrium properties, but will have an 
even more important impact on the evaluation of out-of-equilibrium properties, such as in 
the case of simulated AFM pulling experiments. 

Still, much has to be done regarding the puckering properties of carbohydrates. In par- 
ticular, the role of skew conformations — possibly detected in the NMR experiments of 
Snyder and Serianni^- but not significantly present in both the 45a4 and new parameter sets 
— has to be clarified, and different sugars, such as pyranoses of pentose nature, should be 
investigated to further test the force field. Eventually, the puckering properties of furanoses 
should deserve further attention. The GROMOS force field, however, is not the only one which 
has been found to be problematic in reproducing proper puckering properties, and for many 
force fields the investigations of ring conformer populations are scarce, if not missing at all. 
Our hope is that, beside the usefluness related to the specific case of the GROMOS force field 
parametrization, this work could also serve to attract attention to the importance of the 
puckering problem in carbohydrate simulations, and to stimulate further investigations. 
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e (rad) 

FIG. 3: Puckering free energy landscape (already corrected with the results from the umbrella 
sampling phase) of /3-D-Glc using the standard 45a4 parameter set. Upper panel: the {9,(j)) plane 
is shown, with isolines drawn every 10 kJ/mol starting from the absolute minimum. Lower panel: 
projection of the free energy profile onto the (p = plane. Darker colors correspond to lower 
energies. 
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FIG. 4: Contribution to the total puckering free energy of /3-D-Glc (45a4 parameters) from the 
equihbrium sampling, —k^T In p^ii^{9, (p). An immaterial constant has been added, so that the 
correction term is zero at the position of the free energy global minimum. Isolines are drawn every 
2.5 kJ/mol (upper panel). The projection onto the (j) = plane of the correction term (lower panel). 
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FIG. 5: Inverted chair free energy difference for the (3 (upper panel) and a (lower panel) series of 
aldopyranoses. The three curves refer to the simulation results obtained using the 45a4 parameter 
set (squares) and the predictions of Ref. 



62l (triangles) and Ref. 



63l (circles). Lines are a guide to 



the eye, and error bars are always smaller than the symbols (0.5 kJ/mol). 
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FIG. 6: Results of the unbiased molecular dynamics simulation of /3-D-Glc. Histogram of the popu- 
lation along the 6 puckering coordinate (upper panel) and time evolution of the number of inverted 
chairs (lower panel). 
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FIG. 7: Inverted chair free energy difference for tlie (3 (upper panel) and a (lower panel) series 
of aldopyranoses. The three curves refer to the simulation results using the new parameter set 
(squares) and the predictions of Ref. 



621 (triangles) and Ref. 



63l (circles). Lines are a guide to the 



eye, and error bars are always smaller than the symbols (0.5 kJ/mol). 
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FIG. 8: Puckering free energy landscapes (in kJ/mol) obtained by the combined metadynam- 
ics - umbrella sampling, for the series of /3-D-aldopyranoses simulated using the new force field 
parameters. Isolines are drawn every 10 kJ/mol, starting from the global minimum. 
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FIG. 9: Puckering free energy landscapes (in kJ/mol) obtained by the combined metadynam- 
ics - umbrella sampling, for the series of a-D-aldopyranoses simulated using the new force field 
parameters. Isolines are drawn every 10 kJ/mol, starting from the global minimum. 
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Stereoisomer 


CI 


C2 


C3 


C4 


C5 


;3-D-Glc 


eq 


eq 


eq 


eq 


eq 


;9-D-Gal 


eq 


cq 


cq 


ax 


eq 


/3-D-Man 


eq 


ax 


eq 


eq 


eq 


/3-D- All 


eq 


eq 


ax 


eq 


eq 


/3-D-Tal 


eq 


ax 


eq 


ax 


eq 


;3-D-Gul 


eq 


eq 


ax 


ax 


eq 


^-D-Alt 


eq 


ax 


ax 


eq 


eq 


/3-D-Ido 


eq 


ax 


ax 


ax 


eq 



TABLE I: Orientation of ring substituents (ax: axial; eq: equatorial) for different stereoisomers 
(referred to the ^Ci conformer). 
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Metnod 


AG[a\ 


AG[pJ 


Molecular mechanical^ 


-0.08 


3.77 


Semiempirical^ 


-2.09 


5.43 


NMR§ 


-4.0 


3.4 



" Free energy of t 



Ref. 



62 ; § Ref. 



67 



le a— anomer (kJ/mol); ^ Free energy of the /3— anomer (kJ/mol); ^ Ref. 63 ; if 
A positive value of the difference denotes the ^Ci as the most stable conformer 



TABLE II: Theoretical and experimental estimates of idose inverted chair free energy. 
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Isomer 


AG [^04]" 


Next** 


AG [Next]'' 




AG[TS]'' 


P[^Ci]^ 


P[^C4]^ 


P[N 


ext]'* 




;3-D-Glc 


10.0 ±0.2 


30b 


16.7±0.2 


(1.06,0.10) 


25.9 ±0.2 


98.42(1) 


1.52(1) 


4.6( 


1) 


xlO" 


-2 


/3-D-Gal 


4.6 ±0.2 




21.5±0.2 


(1.11,3.94) 


38.6 ±0.4 


91.2(1) 


8.8(1) 


6.2( 


1) 


xlO- 


-3 


/3-D-Man 


3.0 ±0.2 




23.2 ±0.2 


(1.11,2.44) 


41.2 ±0.2 


80.6(1) 


19.4(1) 


4.3( 


1) 


xlO" 


-3 


;S-D-A11 


21.2±0.2 




27.9 ±0.2 


(1.17,3.40) 


36.9 ±0.3 


99.98(1) 


0.0205(2) 


4.6( 


1) 


xlO" 


-4 


/3-D-Tal 


-2.8 ±0.2 




36.8 ±0.2 


(1.11,3.62) 


47.6 ±0.2 


58.7(2) 


41.3(1) 


3.8( 


1) 


xlO- 


-5 


/3-D-Gul 


14.8 ±0.3 




33.9 ±0.1 


(1.22,3.30) 


43.2 ±0.2 


99.92(1) 


0.08(1) 


5.7( 


1) 


xlO" 


c 
-0 


/?-D-Alt 


15.5 ±0.2 




33.7±0.1 


(1.11,3.94) 


42.4 ±0.2 


99.87(1) 


0.129(2) 


9.4( 


1) 


xlO' 


c 



/3-D-Ido 


13.0 ±0.2 


B25 


35.0 ±0.3 


(1.06,4.47) 


44.9 ±0.1 


99.59(1) 


0.412(4) 


6.1( 


1) 


xlO" 


-u 


a-D-Glc 


3.7±0.2 




25.7±0.2 


(1.11,0.21) 


35.2 ±0.2 


90.85(6) 


9.15(7) 


1.7( 


1) 


xlO" 


-3 


a-D-Gal 


2.1±0.2 




33.0 ±0.2 


(1.11,3.40) 


47.1 ±0.2 


74.9(2) 


25.1(2) 


5.8( 


1) 


xlO" 


-5 


a-D-Man 


2.3 ±0.3 


OS2 


15.2 ±0.3 


(1.11,3.30) 


29.5 ±0.4 


71.9(3) 


28.0(2) 


3.9( 


1) 


xlO- 


-2 


a-D-All 


18.2 ±0.2 


30b 


35.6 ±0.2 


(1.65,4.79) 


56.6 ±0.3 


99.94(1) 


0.060(1) 


3.4( 


1) 


xlO" 


-5 


a-D-Tal 


1.2 ±0.3 




32.1 ±0.3 


(1.11,3.51) 


37.1 ±0.3 


63.7(4) 


36.3(3) 


3.0( 


1) 


xlO" 


-4 


a-D-Gul 


12.9 ±0.2 


OS2 


42.7±0.2 


(1.06,3.72) 


30.5 ±0.2 


99.49(1) 


0.51(1) 


7.0( 


1) 


xlO- 


-6 


a-D-Alt 


18.5 ±0.4 


OS2 


23.1 ±0.2 


(1.17,2.87) 


45.5 ±0.3 


99.93(1) 


0.065(1) 


5.5( 


1) 


xlO" 


-3 


a-D-Ido 


14.0 ±0.2 




27.9 ±0.2 


(1.11,4.36) 


33.6 ±0.2 


99.47(1) 


0.524(3) 


1.4( 


1) 


xlO" 


-3 



Energies are expressed in kJ/mol, and probabilities in percent; " Free energy of ^€4; * Next most 



populated conformer after ^Ci and ^€4; ^ Free energy of the next most populated conformer; Location of 
the transition state {9,(j)); ^ Free energy of the transition state ; Population of ^Ci (%); ^ Population of 
^€4 (%); Population of the next most populated conformer. When a basin is present, which encompasses 
many states, all conformers involved are listed. 

TABLE III: Free energy and population of different conformers using the 45a4 parameter set. 



45 



dihedral 


angle 




cos 5 


n 


C3-C2- 


-Cl-01 t 


0.5 


-1 


2 


C4-C3- 


-C2-02 t 


0.5 


-1 


2 


C1-C2- 


-C3-03 t 


2.4 


-1 


2 


C5-C4- 


-C3-03 t 


2.4 


-1 


2 


Cl-05- 


-C5-C6 * 


0.5 


-1 


2 



^ in kJ/mol; 'f interaction modified with respect to the correspontent 45a4 one; ^ new interaction term (not 
present in 45a4); The functional form for the torsional interaction is that of EqfTOl All other interactions 
are the same as in the 45a4 se1j2^. 



TABLE IV: New parameters for D-aldopyranoses torsional interactions. 
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Isomer 


AG [^€4]" 


Next 


AG [Next] 


rpgd 


AG[TS]^ 




P[lC4]^ 


P[N< 


BXt]'* 




;3-D-Glc 


27.1 ±0.3 


30b 


30.9 ±0.2 


(1.17,3.19) 


45.9 ±0.5 


99.99(1) 


0.002(1) 4.17(4) X 10" 


-4 


;9-D-Gal 


19.3 ±0.2 




35.6 ±0.2 


(1.17,3.51) 


39.7±0.2 


99.95(1) 


0.051(1) 


8.8(1) 


xlO- 


-5 


/3-D-Man 


13.1 ±0.4 




32.9 ±0.4 


(1.17,2.66) 


42.7±0.5 


99.58(1) 


0.42(1) 


1.41(3) X 10^ 


-4 


;S-D-A11 


12.5 ±0.2 




21.8±0.2 


(1.22,3.19) 


36.8 ±0.3 


99.31(1) 


0.68(1) 


9.23(1) 


xlO" 


-3 


/3-D-Tal 


11.9±0.3 


B3O 


43.2 ±0.1 


(1.11,3.40) 


46.8 ±0.2 


99.53(1) 


0.47(1) 


3.0(1) 


xlO- 


-6 


/3-D-Gul 


13.3 ±0.3 


°S2 


31.5 ±0.4 


(1.17,3.19) 


45.9 ±0.5 


99.00(2) 


0.99(1) 


3.4(1) 


xlO" 


A 
-4t 


/?-D-Alt 


6.1±0.1 




24.7 ±0.1 


(1.17,3.72) 


43.8±0.2 


93.88(4) 


6.12(4) 


2.4(1) 


xlO" 


Q 


/3-D-Ido 


4.1 ±0.3 


B25 


27.4 ±0.2 


(1.06,5.85) 


44.9 ±0.2 


88.51(9) 


11.5(2) 


8.6(1) 


xlO" 


A 


a-D-Glc 


16.2 ±0.2 


^Sa-^^B 


35.1 ±0.2 


(1.22,2.66) 


42.5 ±0.2 


99.88(1) 


0.124(1) 


8.6(1) 


xlO" 


-5 


a-D-Gal 


14.5 ±0.2 




44.7 ±0.2 


(1.22,3.19) 


48.5 ±0.4 


99.58(1) 


0.416(5) 4.0(1) 


xlO- 


-6 


a-D-Man 


6.5 ±0.2 


°S2 


22.8 ±0.1 


(1.11,3.19) 


32.3 ±0.3 


95.72(3) 


4.27(5) 


6.0(1) 


xlO" 


-3 


a-D-All 


8.1 ±0.2 


Bso-^Ss 


27.1 ±0.2 


(1.17,3.40) 


38.9 ±0.2 


97.92(2) 


2.08(3) 


8.7(1) 


xlO" 


-4 


a-D-Tal 


6.8±0.1 


0S2 


33.7 ±0.2 


(1.65,4.47) 


52.5 ±0.2 


91.60(6) 


8.40(6) 


1.7(1) 


xlO" 


-4 


a-D-Gul 


4.2 ±0.2 


0S2 


30.3 ±0.2 


(1.49,3.94) 


49.2 ±0.3 


90.44(7) 


9.6(1) 


1.2(1) 


xlO- 


-4 


a-D-Alt 


2.8 ±0.3 


B30 


15.0 ±0.2 


(1.06,3.72) 


28.9 ±0.2 


79.2(1) 


20.7(2) 


8.74(1) X 10- 


-2 


a-D-Ido 


1.1±0.2 


°S2-^''B 20.6 ±0.2 


(1.11,5.53) 


31.6±0.2 


51.8(3) 


48.2(2) 


1.38(1) xlO- 


-2 



Energies are expressed in kJ/mol, and probabilities in percent; " Free energy of ^€4; * Next most 



populated conformer after ^Ci and ^€4; ^ Free energy of the next most populated conformer; Location of 
the transition state {9,(j)); ^ Free energy of the transition state ; Population of ^Ci (%); ^ Population of 
^€4 (%); Population of the next most populated conformer. When a basin is present, which encompasses 
many states, all conformers involved are listed. 

TABLE V: Free energy and population of different conformers using the new parameter set. 
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